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ABSTRACT 

We study time evolutions of superfluid neutron stars, focussing on the nature of 
the oscillation spectrum, the effect of mutual friction force on the oscillations and the 
hydrodynamical spin- up phase of pulsar glitches. We linearise the dynamical equations 
of a Newtonian two-fluid model for rapidly rotating backgrounds. In the axisymmetric 
equilibrium configurations, the two fluid components corotate and are in /3-equilibrium. 
We use analytical equations of state that generate stratified and non-stratified stel- 
lar models, which enable us to study the coupling between the dynamical degrees of 
freedom of the system. By means of time evolutions of the linearised dynamical equa- 
tions, we determine the spectrum of axisymmetric and non-axisymmetric oscillation 
modes, accounting for the contribution of the gravitational potential perturbations, 
i.e. without adopting the Cowling approximation. We study the mutual friction damp- 
ing of the superfluid oscillations and consider the effects of the non-dissipative part 
of the mutual friction force on the mode frequencies. We also provide technical de- 
tails and relevant tests for the hydrody namical model of pulsar glitches discussed by 
Siderv. Passamonti fc Andersson ( 2010[ ). In particular, we describe the method used 
to generate the initial data that mimic the pre-glitch state, and derive the equations 
that are used to extract the gravitational-wave signal. 

Key words: methods: numerical ~ stars: neutron - stars: oscillation - star:rotation 
- gravitational waves 



1 INTRODUCTION 

Mature neutron stars are expected to have superfluid and superconducting components in their interior. Shortly after a 
neutron star's birth the temperature decreases below T ~ 10^ K, at which point superfluid neutrons should be present 
both in the inner crust and the outer core, while the core protons should form a superconductor. At all relevant temper- 
atures, the electrons form a "normal" fluid that is tightly locked to the protons due to the electromagnetic interaction. 
This suggests th at the dynamics of mature neutron stars depe nds on the detailed interaction between coupled superfluids- 
superconductors ( Glampedakis. Andersson fc SamuelssorJ 2O10l ). i.e. represents a complex physics problem. The situation is 
not expected to simplify if one also accounts for the inner neutron star core, at several times the nuclear saturation density, 
where exotic states like hyperon superfluid mixtures or deconfined quark condensates may be present. 

Although it is generally appreciated that neutron stars have this very complicated structure, the evidence for the presence 
of the different superfluid phases remain indirect. The strongest support comes from observed pulsar glitches, rapid spin-up 
events seen in a number of young pulsars (and also some magnetars) during their magnetic slow-down phase. The typical glitch 
size is very small, representing a relative change (Af2) in the observed rotation rate (Q) in the range 10~^ < ASl/Sl < 10~^. The 
currently accepted model for these events relies on the transfer of angular momentum between a (faster spinning) superfluid 
neutron component and the star's (slower spinning) elastic crust (to which the magnetic field is anchored). The exchange 



* E-mail:a. passamonti@soton.ac.uk 



2 A. Passamonti & N. Andersson 



is thought to be mediated by neu tron vortices (by means of which the superfluid mimics bulk rotation) and the associated 
mutual friction (jAlpar et al.lll984f ). 

A challenge for future observations is to probe the detailed physics of a neutron star's interior. In this context, as- 
teroseismology associated with either gravitational or electromagnetic signals seems particularly promising. In fact, the 
quasiperiodic oscillations seen in the tails of giant magnetar flares m ay have provided us with the first opportunity to 
test our theoretical models against observational data (see for instance I Watts fc Strohmaveii 120071 . and references therein). 
The observed variability likely originates from crustal oscillations and depends on the detailed crust dynamics and the 
interaction with the neutron star's magnetic field. These observations have led to a resurgence of interest in neutron- 
star seismology and a renewe d assault on the problem of magn etic star oscillations, a seriously challenging problem from 
the theory point-of-view (see IColaiuda. Beyer fc Kokkota 3 I2OO9I , for a discussion of the literature). In the context of the 
prese nt paper, the potential relevance of the neutron superfluid that penetrates th e neutron star crust is particularly rele- 
vant ( Andersson. Glampedakis fc Samuelsson 20091 : Samuelsson fc Andersson 20091 ). The prospect of detecting gravitational 
waves from oscillating neutron stars is also exciting, especially sin ce the associated signals will allow us to probe the high- 
density region and hence the supranuclear equation of state (EoS) ( Andersson fc KokkotasI 19981 : Benhar. Ferrari fc Gualtieri 
2OO4I : ISamuelsson fc Anderssonll2007l : lAndersson et allEoogl ). 



In order to faciliate future observations and the decoding of collected data, we need to improve our models considerably. 
The superfluid aspects are particularly interesting in this respect, since the oscillation spectrum of a superfluid star is more 
complex than that of a single fluid model. In superfluid regions fluid elements can execute both co- and counter-moving 
motion, leading to the existence of unique "superfluid" oscillation modes. Our understanding of the nature of the additional 
degree(s) of freedom and the effect on observables must be improved by detailed modelling, ultimately in the context of 
general relativistic multi-fluid dynamics. 

The present work presents recent progress towards this goal. We study the oscillations of superfluid neutron stars by 



evolvin g in time the linearized two-fluid equations in Newtonian gravity. We improve on the analysis of iPassamonti et al 



( 2009al l by including the perturbations of the gravitational potential. We also account for the mutual friction force associated 
with vortices, and implement quadrupole extraction of the gravitational-wave signal associated with the fluid motion. We 
provide the detail ed analysis (and rel evant code tests) for the configurations that we recently used to study the hydrodynamics 
of pulsar glitches ( Sidery et al. 2010| ). We consider two simple analytical EoS and construct two distinct sequences of rapidly 
rotating stars, the main difference being the presence or absence of composition gradients. Such gradients impact on the 
superfluid dynamics, as the co- and counter-moving degrees of freedom are coupled in stratified models. From time-evolutions 
of the relevant perturbation equations, with the gravitational potential perturbation included, we determine the axi- and non- 
axisymmetric oscillation modes for models that rotate up to the mass shedding limit. Finally, we account for the (standard 
form of the) mutual friction force. This adds two coupling terms to the equations of motion. One component is dissipative and 
damps an oscillation mode, while the other modifies the frequencies of the superfluid modes. We study both these effects and 
infer an analytical relation for the associated frequency change of the non-axisymmetric superfluid fundamental and inertial 
modes. 



2 EQUATIONS OF MOTION 



In a basic model for superfluid neutron stars, the matter constituents are superfluid neutrons, superconducting protons and 
normal electrons. Given the typical dynamical timescale of stellar oscillations, one would expect the charged particles to be 
efficiently locked together by the electromagnetic interaction. Therefore, the dynamics of superfluid stars depends on two 
components, a neutron superfluid and a neutral conglomerate of protons and electrons. For simplicity, we will refer to the 
lat ter mixture as th e " protons" in the following. More detail e d disc ussion and justification for the two-fluid model is provided 
bv lMendelll (|l991al lbl). |Prb3 (|2004l ) and I Andersson fc Comej (|2006l ). 

When the mass of each fluid component is conserved, i.e. when we neglect the various particle reactions, the dynamics 
of a superfluid star is described b y two mass conservation laws, two Euler-type equations and the Poisson equation for the 
gravitational potential (|Prixll2004l ). These take the form; 



9tpx + Vi ( /5x«x ) = 



{dt + Wx Vfc) (uf + Exwr) + V, ($ Ax) + e^wl^V, 



Px 



(1) 



(2) 



V^-I- = 47rG'/9. (3) 

These equations are given in a coordinate basis, which means that the indices i and k denote spatial components of the various 
vectors. Meanwhile the indices x and y label the two fluid components. In the present case these constituent indices will be n 
for the neutrons and p for the protons. Throughout this work, the summation rule for repeated indices applies only for spatial 
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indices. In equations (H])-©, the total mass density is p = pn + Pp, /Xx is the chemical potential for each fluid component 
(scaled with the particle mass m — rrin — rUp), <1> is the gravitational potential, while the relative velocity between the two 
fluids is w'^^ = vf — . The parameter accounts for the non-dissipative entrainment effect. In a neutron star core the 
entrainment is due to the strong interaction between the nucleons. From equation ((2]|, it is clear that it leads to a momentum 
that is not longer aligned with the individual component velocity. The vector field f represents the force density acting on 
the X fluid component. In this paper we consider only the vortex mediated mutual friction force. The general form of this 
force is 

fi = 2pn (j3' eijk Q.^w'^y + B e,jfc e'lrn O,'' w^y^ , (4) 

where H' = 51'/^ represents the bulk rotation (later we will assume that the two fluids co-rotate in the unperturbed back- 
ground), and B and B' are the mutual friction parameters. 



2.1 Equation of State 



The equation of state (EoS), that is needed to close the system of equations, can be described by an energy functional 

£ = £ (pn,Pp,Wnp) , (5) 

that ensures Galilean invariance. The chemical potential px and the entrainment parameter are then defined by 



Ex 



9px 

2px 



d£ 



dv 



(6) 
(7) 



When the relative velocity between the two fluids is small, as is the case in most systems of practical relevance, equation (O 
can be expanded in a series: 

£ ^ £o (pn, Pp) + ao (pn, Pp) wip + O [wip) , (8) 

This has the advantage that the bulk EoS £o and the entrainment parameter qq can be independently specified at Wnp — 0. 
From equation ([7]) it follows that the entrainment parameter £x is related to the function ao by 

PxEx = 2ao . (9) 

Despite recent developments ( Chamel 2008h . we do not yet have a realistic EoS that consistently describes the superfluid 
properties of a neutron star core. Therefore, we consider two analytical EoS, based on generalisation of the familiar n — 1 
polytrope. These models are particularly useful if we want to explore the role of entrainment, composition stratiflcation 
and symmetry energy. Moreover, since the two EoS have been used elsewhere we have "independent" tests of our numerical 
results. The main difference between our two sets of models is the presence, or absence, of composition gradients. This is 
important since the co- and counter-moving degrees of freedom are coupled in stratifled neutron stars, which means that the 
gravitational- wave spectrum may contain the imprints of "superfluid" modes (see Sec. [Sjl . This would not be the case in a 

non-stratified model. 

The first EoS is determined by the following expression (jPrix et al.ll2002l ; lYoshida fc Eriguchil |2004| ; iPassamonti et al 



2009aD: 



£o 



K 



2Ka 



l-(l-ha)xp" 1- (H-a)xp"''"^P ' Xp[l-{l + a)xp] 
where _ R" is a polytropi c constant, Xp is the proton fraction and a is a pa rameter that can be rela ted to the symmetry 



K[l + a 



il + 2a)xp] 2 

rp 1 



(10) 



energy IjPrfac et al.ll2002l ). In this EoS, both Xp and a are taken to be constant (|Passamonti et al.ll2009al ). Using equation (|10|) . 
we construct a sequence o f co-rotating axis ymmetric conflgurations without composition gradients. These correspond to the 
A models used by Pass amonti et al.l (|2009al ) . 

In order to study the effects of stratiflcation on the oscillation spectrum we consider a second EoS, defined by ( Prix fc RieutordI 
2OO2I : lAndersson et ahlbood : IPassamonti et al.ll2009al ): 

£o = kpZ" + kpp;^ . (11) 

Here, the coefficients fcx and 7x are constants. We consider 7n — 1.9, 7p = 1.7 for all the rotating models. In the numerical 
code, the coefficients fcx are given in units of GR^qPo~^'' , where G is the gravitational constant and Req is the equatorial 
radius of the stellar model. We t ake them to have the va lues kn = 0.682 and fcp — 3.419 for the non-rotating model, which 
corresponds to model III used by Prix Sz RieutordI ( 2002[ ). Note that for rotating models, the dimensionless fcx can assume a 
different value with respect to the non-rotating star. For instance, when we impose that the central proton fraction is constant 
for all the sequence of rotating models (see Section [3. If) . From equations (|6]) and (|lip it follows that the chemical potential 
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Table 1. This table provides the main parameters for the two sequences of rotating models. The first column labels each model. In 
the second and third columns we give, respectively, the ratio of polar to equatorial axes and the angular velocity of the star. In the 
fourth column, the rotation rate is compared to the Kepler velocity Q,x that represents the mass shedding limit. The ratio between 
the rotational kinetic energy and gravitational potential energy T/|iy| and the stellar mass are given in the fifth and sixth columns, 
respectively. Finally, the seventh column gives the value of the chemical potential at the centre of the star. All quantities are given in 
dimensionless units, where G is the gravitational constant, po represents the central mass density and R^q is the equatorial radius. 



Model 


Rp 1 Req 






T/\W\ X 102 


M/ipoRlq) 


fio/iGpoRlq) 


A n 


1 nnnnn 
i.UUUUU 


n nnnnn 
U.UUUUU 


n nnnnn 
U.UUUUU 


n nnnnn 
U.UUUUU 


1 OT'iO 
L.Z 1 OZ 




Al 


0.99792 


0.05913 


0.08153 


0.05802 


1.2701 


1.2701 


A2 


0.98333 


0.16675 


0.22992 


0.38482 


1.2479 


1.2477 


A3 


0.95000 


0.28729 


0.39613 


1.16918 


1.1967 


1.1962 


A4 


0.90000 


0.40268 


0.55524 


2.38295 


1.1186 


1.1179 


A5 


0.80000 


0.55626 


0.76700 


4.93320 


0.9557 


0.9576 


A6 


0.70000 


0.65789 


0.90713 


7.56798 


0.7801 


0.7917 


A7 


0.60000 


0.71733 


0.98909 


9.86465 


0.5794 


0.6176 


A8 


0.55625 


0.72524 


1.00000 


10.2760 


0.4749 


0.5361 


CO 


1.00000 


0.00000 


0.00000 


0.00000 


1.0826 


1.1755 


CI 


0.99792 


0.55856 


0.08403 


0.04561 


1.0798 


1.1725 


C2 


0.98333 


0.15764 


0.23716 


0.36682 


1.0601 


1.1516 


C3 


0.95000 


0.27145 


0.40837 


1.11303 


1.0146 


1.1034 


C4 


0.90000 


0.38006 


0.57177 


2.26285 


0.9447 


1.0301 


C5 


0.80000 


0.52334 


0.78732 


4.64841 


0.7984 


0.8792 


C6 


0.70000 


0.61536 


0.92576 


7.01965 


0.6392 


0.7223 


C7 


0.60000 


0.66249 


0.99667 


8.77258 


0.4562 


0.5559 


C8 


0.57656 


0.66471 


1.00000 


8.87526 


0.4077 


0.5147 



and mass density are related by 



fcx7x 



(12) 



where the polytropic index is given by A^'x = (7x — 1) ^ ■ From this result we can determine the proton fraction for a given 
stellar model by imposing /3-equilibrium. After some calculations, we obtain: 



1 



(7pfcp) 



(7nfcr 



--'Vn 



(13) 



This EoS will be used to construct a sequence of stratified rotating models, where the central proton fraction is fixed, 
^p(O) = O-l- These equilibrium configurations have already been used by Siderv et al. ( 201Cll ). and will be refered to as models 
C in the following. 



3 EQUILIBRIUM CONFIGURATIONS 

We study the oscillations of rotating axisymmetric background models where neutrons and protons are in /3-equilibirum and 
co-rotate with constant angular velocity, i.e. we have fin = fip- In this work, we also assume that two fluid components 
coexist throughout the stellar volume. This is obviously artifical; the outer region of a real neutron star will not be superfiuid. 
However, at this stage our main interest is in the bulk core dynamics. In the future we plan to extend our model to account 
appropriately for the expected superfluid regions. At that point we will also consider the role of the elastic crust. 

In Sec. 13.11 we introduce the equations that govern stationary co-rotating equilibrium configurati ons and solve them 



numer ically for the EoS (|10|) and (|11|) . In Sec. 13. 21 we then describe the perturbative approach developed bv lYoshida fc Eriguchi 



1 20041 ) for determining stationary configurations in which the two fluids rotate with a small velocity lag. With this method 
we can obtain non-corotating models as small deviations from a co-rotating equilibrium. Subsequently, we use this approach 
to determine the initial conditions for hydrodynamical glitch evolutions. 



3.1 Corotating background 

The equations that describe rapidly and uniformly rotating background models can be derived by im posing the condi- 
tions of stationarity and axi-symmetry on the Euler-type equations ((2]) and the Poisson equation ((3| I Prix et al. 20021 ; 
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Yoshida fc Eriguchill2004 ). This leads to 



$(r) 



C'x ; 

P(r') 



■r' 



-dr' . 



(14) 
(15) 



where and Cx are, respectively, the angular velocities and the integration constants for the neutron and proton fluids. For 
corotating background models, i.e. when — flp = i}, in which the two fluids are in /3-equilibrium and share a common 
surface, the hydrostatic equilibrium equation p4[) becomes 

„2 



r ■ 2 

M + $ - y sm ( 



n=c, 



(16) 



where fl = jlp — fin is the background chemical potential and C = Cn = Cp. When the system of equations (|15|l - (|16p is 
closed by an EoS, we can numerically determine a corotating stationary axisymmetric background via the self-consistent held 
method (jHachisulllQS^ : iPassamonti et al.ll2009al). The solution is suc h that the surface of the star corresponds to the zero 
chemical potential surface, jl {R (6) ,9) — (jYoshida fc Eriguchill2004l ). 

For the EoS (|10[) and (|11[) we construct the sequences of rotating models A and C, respectively. The set of models extends 
from a non-rotating model up to the mass shedding limit. In the numerical code, we re- write the background equation in 
dimensionless form by using the gravitational constant G, the central mass density po and the equatorial radius Req- All stellar 
models of the sequence have the central proton fraction set to a:;p(0) = 0.1. By specifying the axis ratio between the polar and 
equatorial radius Rp/R^q, the iterative numerical routine determines all other quantities of an axisymmetric configuration. 
The main properties of the rotating models are given in Table [T] From these quantities we can easily construct stellar models 
in physical units. For instance for models C, we can evaluate equation H12|) at the centre and obtain: 

7nfcn = (1 - a;p(0))^~^" , 7pfcp = /i;xp(0)^"^P , (17) 

where the asterisk denotes the dimensionless quantities k* = k^/ {GR^qpQ~^") and fl* = fij (^GpoReq)- Combining equa- 
tions H17p with the dimensionless mass M* = M/ [poReq), we can derive the equatorial radius of the star. 



Req 



kn-Jn 



G 



Ar"-'(l-xp(0))' 



n l/(37„-4) 



(M 



*7n — 2 



-l/(37n-4) 



(18) 



where the physical mass M and the EoS parameters can be arbitrarily chosen. The central mass density, po, and the rotational 
period, P, are determined by the following equations: 



po = 2.8 X 10''' (Af )" 

P = 0.4596 (A//*)^''^ {Q.* 
where £7* = Q./y/Gpo. 



M 



1AM, 



( Re 



M 



\10 km 

-1/2 



1AM, 



g cm 



Req 

10 km 



3/2 



ms , 



(19) 
(20) 



3.2 Non-corotating solutions 



In a multi-fluid system, like an astrophysical neutron star, the various fluid components can have different velocities. This 
is, in fact, an essential element in the favoured model for pulsar glitches where the sudden observed spin-up is explained 
as a transfer of angular momentum between an interior superfluid neutrons and the charged component. In this model, the 
momentum transfer is due to the interaction between the crust and an array of quantised neutron vortices that are generated 
by the stellar rotation. During the magnetically driven spin-down of a neutron star, these vortices are pinned to the crust and 
corotate with the charged components. Therefore, a velocity lag develops between superfluid neutrons and the crust and an 
increasing Magnus force acts on the vortices. When this force becomes stronger than the pinning force, the vortices should 
unpin. At this point they are free to move and can accelerate the crust, generating a glitch. 

Typically, the spin variation observed in a glitch is very small, 10"^ < AQ/n < 10"^ This means that the effects 
of a glitch on the stel l ar st ructure is expected to be tiny and can be studied perturbatively. The approach developed 
bv lYoshida fc Eriguchil (|2004l ) is particularly appropriate for this kind of problem, as the non-corotating quantities are con- 
sidered as small deviations from a stationary, rapidly corotating configuration. We have already us ed these non-corot ating 
corrections as initial data for studying the post-glitch dynamics and the associated stellar oscillations I Siderv et al. 2O10l ). We 
will now provide f urther details about the m ethod. 

Adopting the lYoshida fc Eriguchil (|2004l ) approach, we expand equations (|14|l - (|15p up to the first order in 



(a.-f7p)/(if7nl + lf7p 
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Figure 1. In this figure, we compare our numerical results to the analytical solution of lPrix et al.l ||2002| ) for two slowly rotating stellar 
models with axis ratio 0.996 (left panel) and 0.975 (right panel). These two background stars are described by the EoS l llOll and have 
the same proton fraction Xp = 0.1 and symmetry energy term cr = 0.5. The non-corotating corrections are determined by choosing the 
relative angular velocity (<5f2n , <5f2p ) = (1,0) and imposing the constant central chemical potential condition H28|l . In the two panels, we 
show the radial profile of the perturbed neutron mass density 5pn/po for three different angular directions, i.e. 6 = 0, 7r/4 and 7r/2. Our 
numerical results (solid line) agree very well with the analytical solution (empty circle) for the slowest rotating model (left panel). For 
faster rotating models the slow-rotation solution is expected to be less accurate. This is already evident for the case in the right panel, 
where the numerical and analytical solutions start to disagree. 



Thus, we have 

Px — Pc + (5px , 



(21) 
(22) 
(23) 
(24) 



where the subscript "c" denotes the corotating values. Note that by definition 50.^ represents the relative deviation of the x 
fluid angular velocity with respect to the corotating background, i.e. 50,^ — (Six — f2c) /i^c- For the non-corotating corrections. 
Equations (fT4l) - (fT5|) become 



5px + 5$ - sin^ 9 O.^ Sn^ = 5Cx 



5$ (r) = -G 



- 5p(r') 



dr'. 



The system of equations is closed by 



X- 9px 
Opx = — 



Pn Pp 



Spy, , 



(25) 
(26) 

(27) 



that relates the mass density and the chemical potential perturbations for co-rotating backgrounds. 

Non-corotating solutions can be constructed with either fixed central chemical potential or total mass. For the first class 
of models, we can impose the condition (5/ix|^^o = at the star's centre, and determine the integration constant 5Cx from 
equation ([251 (jVoshida fc Eriguchilliooi ): 



(5Cn — SCa 



For solutions with constant mass, we impose a constraint on the mass of each fiuid component, i.e. 



5Mx 



dr Spx — . 



(28) 



(29) 



In equation (|25|) . we can replace the chemical potential by the mass density perturbation using equation (|27|) . and integrate 
over the star's volume V. The integration constant 5Cx is then given by the following expression: 



<5Cx 



drpx"^" = J drpx"^" (5$ - sin^ enjD.^ 



(30) 
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For the EoS (|10p . the adiabatic index is 7x = 2 and the boundary condition l|30p therefore reduces to: 

SC^^ ^ J dr {S^ - sin^ enjQa:) ■ (31) 



The system of equations (|25p - (|26|l can be solved iteratively. First of all, for a given EoS we determine the co-rotating 
background with the self-consistent field method of Hachisu (198(i), where we specify the axis ratio of the star. Secondly, we 
choose the relative angular velocity of each fluid component. The iteration algorithm then proceeds as follows: i) we solve 
the perturbed Poisson equation (|26|) for an initial guess of the perturbed mass density Sp, ii) we get the integration constant 
5Cx imposing either the condition (|28p or pO[) . Hi) we determine the chemical potential Sfi^ from equation (|25|) and then the 
new mass density 5p from the EoS. This procedure is iterated until the difference between the quantities is smaller than a 
prescribed error. 

An important property of this linear perturbation approach is that we can construct two independent solutions to 
equations (|25p - (|26ll . respectively corresponding to {50,-^,50,^) = (1,0) and ((5r2n,5np) = (0, 1). Since the problem is linear, 
any non-corotating configuration can be obtained as a linear combination of these two solutio ns. 

We have tested our code against the analytical solution for the EoS (|10p determined by Prix et al. ( 2002l l in the slow- 
rotation approximation. We select two slowly rotating models with axis ratio 0.996 and 0.975, respectively. The models have 
the same proton fraction and symmetry energy term, i.e. = 0.1 and a = 0.5. The non-corotating corrections correspond 
to a relative angular velocity {5Q,n, SQ,p) — (1,0) with constant central chemical potential, c.f. (|28|) . In Fig. [Jl we show the 
radial profile of the perturbed neutron mass density Spn/po for the three angles 9 — 0,7r/4 and n/2, respectively. In the 
slowest rotating model, the agreement between the numerical and the analytical solutions is evident. In the second model, 
with axis ratio 0.975 the two solutions begin to differ, as expected. The slow-rotation solution becomes less accurate as the 
star's rotation increases. The same behaviour is found for non-corotating solutions with constant mass, i.e. when SM^ — 0. 
This comparison gives us confidence in our numerically generated background models. 

For the sequence of constant mass models, we show in Figs. [2]and[3]the non-corotating mass density pertubations Sp^/po 
and the gravitational potential perturbation J$ for the C2 model. These are solutions to equations (|25p-(l26l) with SQp = — 10~^ 
and 50,n = 7.74 x 10~*, which were used as initial conditions for the glitch simulations discussed by Siderv et al. ( 201Cll l. 



4 PERTURBATION EQUATIONS 

The dynamics of a superfiuid neutron star can be studied by linearizing the system of differential equations (O"©. In the 
inertial frame, the Eulerian perturbation equations are given by 

dt (5vx + SxcSwyx) + !^ dsSv^ = - V5/ix - V5$ -2nxSv^ + , (32) 

{dt + nd^) Sp^ = -V • (px(5vx) , (33) 
V^5<I> = AivGSp, (34) 
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where is the azimuthal angle associated with the rotational motion, and the perturbed mutual friction force is (in the case 
of a co-rotating background) 

5r = 2pnS'r2 X 5wxy + 2pnB S7 X n X 5wxy . (35) 

The chemical potential perturbations can be expressed in terms of the mass density perturbations using equation (|27p . 

In order to solve numerically equations (|32p - p3[l we use the conjugate momentum perturbations iSpx as dynamical 
variables. These are given by 



(5pn = (1 - £n) <5Vn + £n(SVp , 

(5pp = ep(5vn + (1 - ep) (5vp , 

where we recall that PpEp = Pn£p. By inverting these relations we can determine the velocity fields at any time step, 



(5Vn = 



(1 — Ep) (5pii — £ii5pp 

-£p(5pn + (1 - £n) I^Pp 
l-£ 



(36) 

(37) 



(38) 
(39) 



where £ = £n + £p = £n/a;p. 

The time evolution of the non-axisymmetric perturbation equations is a three-dimensional problem in space. However, 
linear perturbations on an axisymmetric b ackground can be expanded in terms of a set of basis functions (cos mcj) , sin m(j)) , 

198C ). The mass density pertu rbations as well as the other 
2OO2I : IPassamonti et aLll2009al ) 



where m is the azimuthal harmonic index i Papaloizou fc Pringlel 
perturbation quantities then take the following form jjones et al 



6p{t,r,6,(p) = [Sp^ {t,r, 9) cos m(l> + Sp^ {t,r, 9) smrrKj)] 



(40) 



With this Fourier expansion the perturbation equations decouple with respect to m and the problem becomes two-dimensional. 
In particular, for the axisymmetric case (m — 0) only the Sp^ component survives. 



4.1 Boundary Conditions 

In this work, we study axisymmetric (m — 0) and non-axisymmetric oscillations (m = 2) of a superfluid neutron star with 
equatorial and rotational axis symmetry. The numerical domain extends over the region ^5 r/R{6) ^ 1 and ^ 9 7r/2, 
and we need to impose boundary conditions at the surface, origin, rotational axis and equator. 

We first discuss the boundary conditions at the origin (r = 0) and the rotational axis (^ = 0), where the perturbation 
equations must be regular. Let us denote by S^ji a general scalar perturbation, such as the mass density 5p^, the chemical 
potential Sfiy^ and the gravitational potential 5$. For axi-symmetric and non-axisymmetric oscillations, we have to impose 
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the following conditions, respectively 
dSiP 



dr 



r=0 
r=0 



5M 



— for m — , 



= for m = 2 . 



(41) 
(42) 



For the velocity fields 5vx, we impose that there must be no mass flux across the origin (r = 0) for both axisymmetric and 
non-axisymmetric perturbations: 



At the rotational axis 
dSv^ 



0), we impose the following conditions 



0. 



= Sv^ — 5v^ = for m = , 
5v^ = 5v^ = Sv^ = for m = 2 . 



(43) 

(44) 
(45) 



At the equator {6 = 7r/2), the reflection symmetry divides the perturbations into two sets with opposite parity (jPassamonti et al 



20090 ). In the Type I parity class, the scalar perturbations 5^ and the velocity satisfy the following conditions: 

dS^P 



a5< g dSvt „ 

— -— — ov^ — ——— — . 



Meanwhile, the Type II class is such that: 



5ip = 5v^ 



de 



= 5< = . 



(46) 



(47) 



The outer layers of a mature neutron star form an elastic crust made up of nuclei. The crust is an important aspect that 
is yet to be implemented in our numerical model (although we are making progress on it). Our current model is simplified, in 
the sense that we assume that superfluid neutrons and protons are present throughout the stellar volume. We then impose 
the standard boundary condition of a free surface, i.e. require that the Lagrangian perturbation of the individual chemical 
potentials vanish at the surface, i.e. 



A/ix = (5/ix + • V/ic = . 



(48) 



The vector field ^x is the Lagrangian displacement of the x-fiuid component (jAndersson. Comer fc Grosartll2004l '). The value 
of the perturbed chemical potential Sjl^ at the surface is determined from equation (|48[) at each time step. 



5 GRAVITATIONAL- WAVE EXTRACTION 

In order to study the gravitational-wave signal emitted by pulsating superfluid neutron stars, we have implemented the 
quadrupole formula for both axisymmetric and non-axisymmetric oscillations. We will now discuss this implementation, in 
particular, the momentum and stress formula that we use to improve the numerical gravitational- wave extraction. 
The gravitational- wave strain can be determined using the quadrupole formula (Thome 1980l ): 



,2m _ \^ '^^-^ 

m — — l 



(49) 



where i,, is 

''J 



the pure spin tensor harmonic which has "electric-type" parity, i.e. (—1)' ( Thornel 198d ). In this work, we 
focus only on the m — and m — 2 pulsations. In the orthonormal basis of spherical coordinates, the components of the 
(l,m) = (2,0) and {l,m) — (2,2) spin tensor harmonics are, respectively, given by 



1 /15 . 2. 

-W — sm e, 

O \ TV 



(50) 



and 



a. 



= T, 



i 

16 



= — W — (1 + cos e)e 
16 ^ 



■ cos 9 e 



2i<t> 



The quantity X^"' is the quadrupole moment, in the case of a two-fluid star defined by; 



167r 
T5" 



(51) 
(52) 

(53) 
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where the spherical harmonics Y2m for the m = and m = 2 cases are given by 



Y20 = 



Y2: 



1 /5 



4 V TT 



1 /5 



2 V TT 



P20 (cose) 



i./i^sin^0e^'^ 
4 V 27r 



(54) 
(55) 



where P20 (cos 6) is the Legendre polynomial. 

It is well-known that, the numerical calculation of the s econd order time derivative of the quadrupole moment in equa- 
tion (|49p could lead to inaccurate results ( Finn fc Evans|[l990l ). However, the accuracy of the gravitational- wave extraction can 
be improved by transforming equation (|49|) into either the perturbed momentum formula, with a first order time derivative, 
or the perturbed stress formula, where the time derivatives are absent (jPinn fc Evan!j|l990l '). In this work, we use both these 
prescriptions in order to check the wave extraction accuracy. 

For axisymmetric oscillations, m = 0, the gravitational strain can be written as follows: 

^" = §^E<' (56) 



where the quantity A^' is defined by 



Ai 



dr Spx r P 



(57) 



We can reduce the order of the time derivative by using the method developed by Finn fc EvansI ( 199Cl ). and obtain the 
perturbed momentum formula: 



Ai" = 2^ / drrp^ I 5vl + 



dt 



2 de 



(58) 



and the perturbed stress formula: 

At" = 2 j drl^-Qr sine p^Svt-^{nr sine fSp^ 



+ — 

4-K 



Vr$Vr(5<I>P (cos e) + Vg^VoS^ P^\sin 6*) + - {Vr^VgS^ + Vg^VrS^) dgP^\cos e) 



(59) 



where the gradient components Vi in eqaution ((59} are determined in the orthonormal spherical basis, i.e. V = (dr, ^dg, -;r^^d(p) ■ 
At the end of the day, the quantity A^" in the strain equation (I56p can be determined from either of the three equations (|57p - 
(ESll. 

For non-axisymmetric oscillations with I = m = 2, the two independent polarizations of the strain can be written as 
follows: 



7 22 -7 22 I 22 ■\r22 

hge - ihg^ = h -2Y , 



where -2Y^^ is the s 



-2 spin- weighted spherical harmonics 

1 /5 



2Y' 



8 V vr 



(l-fcos6')''e''''* 



and we have defined the quantity 



We can then re-write equation (|62[) as follows 



where 



15 r dt2 



drSpr Y2: 



y22 _ ^ ^ ^22 
^ 15 r ^ " ' 



Ai 



df^ 



J dr Sp^r^Y22 ■ 



(60) 



(61) 



(62) 



(63) 



(64) 



In equation (|64p . the order of the time derivatives can be reduced by using the equations of motion (see Appendix |A] for more 
details). This leads to the following expression: 
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t(Gp„) 



t(G p„) 



Figure 4. We compare the gravitational-wave extraction results for axisymmetric and non-axisymmetric oscillations. The signal is 
generated by perturbing the stellar model C2, and the illustrated quantities are dimensionless. The left panel shows the dimensionless 
code quantity determined from three equivalent equations, respectively, the second time derivative of the quadrupole moment (dot- 
dashed line), the momentum-formula (solid-line) and the stress formula (dashed- line). In the right panel, we show the waveform of the 
m = 2 non-axisymmetric oscillations for the perturbed C2 model. The upper and lower right panels displays respectively the real part 
and the imaginary part A-j-^ of the dimensionless quantity A^l^ determined by the code. We compare the signal extraction to the 
momentum-formula (solid-line) and the stress formula (dashed-line) . 



For linear perturbations on a corotating background, we can further transform equation (|65[l into the foUowing expression: 
1 ^15" 



= 2 y 2~ J l"^'""^'" ^i'l ^ ['^'^x + i ^sin ^^Ux -f- cos 61(51)^^1 — (S7r sin S)^ (5/9x 

+ — fsin^ e V,.$Vr5$ + cos^ e Ve$Ve5<l' + sin 9 cos (Vr$Ve(S$ + Ve$V,.(S$) 
- i (sin^Vr* + cosSVe*) V^^*]} , 



(66) 



where the time derivatives are absent. In equations (|65p and (|66p . the perturbations are determined in the inertial frame. 
The energy radiated as gravitational waves is determined by the following equation (Thome 198Cll ): 



p2m _ 1 G 

SZTT & 



dt. 



By using Parseval's Theorem we can write equation (|67p for the (2, 0) and (2, 2) components as follows: 



(67) 



/?22 _ 
^rad 



oo J 771 20 



dEf 

dv 

oo ^^22 



dv ■ 



16 aG 
15 







, 64 3 G 
du 75 C 



roc 

/ 


^2„ 


2 

du , 


Jo 






roc 

/ 


^22 


2 

du , 


Jo 







(68) 
(69) 



where A^"" = A^™ ^ 

The characteristic strain of the gravitational- wave signal is then given by ( Flanagan fc Hughe j 199"^ ): 



An + Ap, and yl is its Fourier transformation. 



he (u) 



2G 1 dE_ 
n'^c? d\du^ 



(70) 



where d is the source distance. The strains fe^" and h^^ are related to the dimensionless quantities Af and Al^ used in the 



numerical code by the following expressions: 



= 1.414 X 10~'''^f(M*)~^ 



= 4.109 X 10""Af (M*)" 



M 



IAMq 
M 



Re 



10 km 



^ / 1 kpc 



1.4Mt 



Re 



10 km 



d 

^ f i kpc 



(71) 
(72) 
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Table 2. Comparison of the first three I = 0,2 ordinary and super fluid mode fr equencies and the lPrix &: Rieutor results. The 

star is the non-rotating CO model, which corresponds to model III of lPrix fc Rieut ord (2002), where the entrainment parameter e is zero. 
Frequencies are given in units of <j/y/Gpo and have been determined with an FFT of the gravitational strain time evolution. For this 
specific numerical simulation the frequency error bar is Ac/ y/GpQ = 5.58 X 10"'^. The Prix and Rieutord values are denoted by PR. In 
the final column we show the relative error between our and PR results. 



Mode a/^/Gpo a/^/GpQ l\cr/cF 
PR [ % ] 



F° 


1.91361 


1.92743 


0.7 


F= 


2.52823 


2.53376 


0.2 


H? 


3.94917 


3.94911 


< 0.1 


Hf 


4.20552 


4.20420 


< 0.1 




5.61069 


5.52870 


1.5 




5.93799 


5.92165 


0.3 


f ° 


1.33511 


1.33178 


0.2 


f 


1.83142 


1.82281 


0.5 


P? 


3.47686 


3.48786 


0.3 


p! 


3.68465 


3.69878 


0.4 


Pi 


5.24187 


5.25802 


0.3 


pl 


5.51946 


5.52876 


0.2 



Similar relations provide the characteristic strain 

3/2 / „ s 1/2 



-3/2 ( M ( Re^V'^ /Ikpc 



= 1.513 X 10"^** 

Kf ^ 2.398 xlO-|^fl(Arr-^^J i^JlukJ' 

As a first test of the numerical implementation, we compare the gravitational-wave extraction formulae for the axisym- 
metric and non-axisymmetric oscillations. We evolve the C2 model with a density perturbation and extract the signal using 
equations H57|) - (I59|I for the m = Q pulsations, and (|65|l - (|66p for the m = 2 oscillations. Typical results are shown in Figure [l] 
We generally find good agreement between the different numerical results, although we note that (as expected) the momentum 
and stress formulae produce a smoother signal than the "raw" quadrupole formula (1571) 



As an additional test, we have used the relativistic numerical code developed by Nag ar fc Diaj ( 20041 ) and Passamonti et al 



( 20071 1 to test the results of the gravitational- wave extraction routine. In the relativistic cas e, the linear perturbations of non- 
rotating relativistic stars were evolved and the signal was extracted using the Zerilli function (|Zerillilll970l ). From the Newtonian 
approach used in the current work, it is evident that we cannot accurately reproduce the relativistic results. However, we 
can establish that our calculations provide a good estimate of the amplitude of the gravitational-wave strain. To this end, 
we consider a star with mass M — IAMq and radius R = 14 km, and evolve the relativistic code with an initial enthalpy 
perturbation, which produces an averaged pulsational kinetic energy of (Ek) — 5.62 x 1O~^M0C^, where c is the speed of 
ligth. The related gravitational-wave strain is almost monochromatic and for a source at lOkpc the maximal amplitude is 
h'^°\ ~ 2.18 X 10"^^ sin^ 

1 max 

With our 2D Newtonian code, we then evolve in time non-radial oscillations of a superfluid non-rotating star for both 
the EoS (|1Q[) and ()lip . The kinetic energy of oscillating superfluid stars can be determined by the following expression: 

Ek = / dr [pn (1 - £n) \5Vn\^ + 2pnen(5Vn ' 5Vp -|- Pp (1 - Ep) |5vpH . (75) 



2 

If we evolve oscillations that have the same pulsational kinetic energy as in the case studied with the relativistic code, we 
obtain h^^l ~ 1.55 x 10"^^ sin^ for model AO and ~ 1.433 x 10"^^ sin^ 9 for model CO. In this calculation, we 

I max I max 

have used equation ()7ip with the parameters of the relativistic stellar model. This test shows that we can be confident that 
the implementation of the quadruple formula in our code provides reasonable results, in accordance with the expected relation 
between the pulsational kinetic energy and the gravitational-wave strain. 



6 RESULTS 

Having formulated the time-evolution problem and described our implementation of the gravitational-wave extraction, we 
will now discuss our results. In this section, we focus on the effects of the gravitational potential perturbation and the 
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Figure 5. This figure displays tiie effect of rotation on the quasi-radial and axisymmetric quadrupole modes. We use the sequence of 
models A, with entrainment parameter e = 0.5, proton fraction Xp =0.1 and vanishing symmetry energy term. On the horizontal axis 
the angular velocity is rescaled with the Kepler angular velocity Qk, while the mode frequencies are given in dimensionless units and 
for a rotating frame. In the left panel, we show some "ordinary modes" , which are due to the co-moving degrees of freedom. We identify 
the / = and I = 2 fundamental modes and the first three quasi-radial overtones and I = 2 pressure modes. In the right panel, we show 
instead the "superfluid modes" , which correspond to the counter-moving degrees of freedom. In this case we show the modes up to the 
second overtones. For these non-stratified models, the ordinary and superfluid modes arc decoupled. 



mutual friction force on axisymmetric and non-axisymmetric oscillations. We also provide a more detailed analysis of the 
gravitational-wave signal generated by the basic glitch model that we discussed in a previous work jSiderv et al.ll2O10l ). 

The pulsation dynamics is studied with a numerical code that evolves in time the system of hyperbolic perturbation 
equations p2l) - (|33|) . solving at each time step the perturbed Poisso n equation (1341). T he part of the code that evolves the 
hyperbolic equations uses the same technology as in previous work ( Passamonti et all 2009a, b). whereas the elliptic equa- 
tion p4l) is solved using a pseudo spectral method. The numerical grid is two-dimensional and covers the volume of the star, 
i.e. the region r ^ Ri^) ^^nd ^ 9 < tt/2. The implementation uses a new radial coordinate x = x{r,6), which is fitted 
to surfaces of constant chemical potential. This allows us to consider stars that are highly deformed by rotation. The pertur- 
bation variables are discretized on this grid and updated in time with a Mac-Cormack algorithm. The numerical simulations 
are stabilised from high frequency noise with the implementation of a fourth order Kreiss-Oliger numerical dissipation. More 
technical details have been discussed in Passamonti et al. ( 2009al lbh . 

In order to solve the elliptic equation (|34p with a spectral method, and save computational time, we set up a second numeri- 
cal grid with lower resolution. This is important, since the spectral solver must be used at each time step, leading to a significant 
slow-down of the simulations. However, the lower resolution on the spectral grid does not affect the r esults, as spectral elliptic 



solve rs provide highly accurate and rapidly convergent solutions already for relatively coarse grids (jGrandclement fc Novak 



Therefore, at each time step we first fit the mass density perturbation Sp on the spectral grid and then use the spectral 



routines to determine the gravitational potential perturbation S^. Subsequently, we fit the new value of 5^ to the original 
grid for the hyperbolic equations and carry on the evolution. The numerical code provides stable simulations for all rotating 
stellar models considered in this paper. 

In this work, our choice of variables differs from that of lPassamonti et al.l (|2009al ). We evolve the velocity perturbations of 
the two-fluids components instead of the "mass flux" perturbations of the co-moving and counter-moving degrees of freedom. 
The two formulations are obviously mathematically equivalent, but we wanted to develop a code based on the new set of 
variables in order to explore which formulation is best suited for future extensions. This is important, as we plan to add more 
realistic physics to our models by implementing an elastic crust reg i on. As a first test, we compare the results of the new code 
to those obtained in Cowling approximation bv iPassamonti et al.l (|2009al ). Neglecting the perturbation of the gravitational 
potential, i.e. setting (5$ = 0, we find a complete agreement between the two numerical codes. 



In order to study the spectral properties discussed below, in Sec. l6.1l and l6.2l we consider "generic" initial conditions that 
excite a large set of oscillation modes. For Type I perturbations we provide the following expression for the mass density: 



Spn — —Spp = 



R{9) 



Yn {6, 0) 



(76) 



where neutrons and protons are initially counter-moving. For Type II perturbations, we excite mainly normal and superfluid 
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Figure 6. In this figure, we show the axisymmetric modes for the sequence of C models determined in a rotating frame. These stellar 
models are stratified and the ordinary and superfluid degrees of freedom are coupled. We identify some of the acoustic modes and their 
dependence on the rotational rate Q/Qk- 



-modes with the foUowing initial data: 



SVn 



-5Vn 



R{e) 



(77) 



where Yjf (6, (j)) is a magnetic spherical harmonic I Thornelll98(t ). For the glitch simulations, we use the non-corotating solutions 
derived in Sec. 13.21 

We test the elliptic solvers by comparing th e mode frequencies extracted from our time evolutions to those obtained in 
the frequency domain by Prix fc RieutordI ( 20021') . We determine the oscillation frequencies of the non-rotating model CO (see 
Table [!}, which corresponds to model III of Iprix fc RieutordI ( 2002 ). For the zero entrainment case, i.e. when (e — 0), the 
results in Table [2] show that the fr equencies determ i ned w ith our code (by an FFT of the time-evolved perturbations) agree 
very well with those calculated bv lPrix fc Rieutord 1I2OO2I) . 



6.1 Spectrum 



The oscillation spectrum of superfluid rotating neutron stars contains the imprints of two-fluid dynamics and of the mutual 
friction force. For a singl e fluid star, th e general mode classification is based on the main restoring force that acts on the 
displaced fluid elements (|Cowlinj|l94ll ). For nonrotating models without magnetic fleld and crust, the spectrum is formed 
by the acoustic, the fundamental and the gravity modes. The acoustic modes are mainly restored by pressure variations and 
cover the high frequency range of the spectrum, above 1 kHz. At lower frequencies, typically below 100 Hz, composition 
and thermal gradients generate the class of gravity modes that are restored by buoyancy. The fundamental mode, whose 
frequency scales with the average stellar density, separates these two classes of modes. In rotating stars, the Coriolis force 
provides an additional restoring force, leading to the presence of inertial modes. Since the frequency of these modes scales 
with the rotation rate, they typically lie in the same low frequency region as the g-modes. For rotating stars with composition 
gradients, the inerti al and gravity modes form a unique class with mixed properties, referred to as gravity-inertial modes (for 
a recent analysis see Passamonti et al. 2009bl : Gaertig fc Kokkotas 20091 ) . 

In addition to this general classification, any oscillation mode can be labeled by the indices {l,m) associated with the 
spherical harmonics Y/"{8,(j)). In spherical stars, this is due to the decomposition of the perturbation functions in vector 
harmonics. For rotating stars, we can use the same description as long as we can track a mode back to its non-rotating limit. 
Finally, for any value of {l,m), the oscillation modes can be ordered by the number of radial nodes in their eigenfunctions. 
The fundamental mode 'f does not have radial nodes, while the series of pressure 'p^ and gravity modes 'g^ have i nodes. 

In superfiuid neutron stars, the additional degree of freedom enriches the dynamics. The two fluids can oscillate both 
in phase and counter-phase. The co-moving degree of freedom produces the class of "ordinary modes", very similar to the 
single fluid resu l ts described above. There i s, however, one import ant difference: the gravity modes are absent in superfluid 
stars l|Lej|l995l : lAndersson fc Comer! I2OOII : iPrix fc Rieutord 120021 ) . The counter-moving degree of freedom generates a new 
class of acoustic and inertial modes, known as "superfluid" modes. These modes strongly depend on the superfluid aspects. 



Hydrodynamics of superfluid neutron stars 15 




such as entrainment and mutual friction. We will label ordinary and superfluid modes by an upper index, for instance the 
I — 2 fundamental ordinary mode will be expressed as ^f °, while ^f represents the corresponding superfluid mode. 

We focus our attention on the quasi-radial {I = 0) and quadrupole (/ = 2) oscillation modes and study their behaviour 
in rapidly rotating models all the way to the mass shedding limit. In non-rotating models, the I = modes are purely radial 
and do not generate gravitational radiation. However, due to coupling of the different multipoles, this is no longer true in the 
rotating case. The quasi-radial fundamental mode will be denoted by F and its i overtones by H^. The quadrupole modes 
(l = 2) are expected to be dominant in the gravitational signal, and we study both axisymmetric and non-axisymmetric 
oscillations. These correspond to m = and m = 2 respectively. 

We start by considering the axisymmetric oscillations for the two sequences of rotating models A and C. For a small 
velocity lag between the two fluids, the entrainment parameter e can be chosen independently fro m the backgro und model 
(see Section |2.1|I . Recent work suggests that it can assume values in the range 0.2 ^ e ^ 0.8 I Chamel 20081 ). Here, we 
consider o nly the case e = 0.5, as the effect of this parameter on the oscillation frequencies has been already discussed 
elsewhere ( Prix fc RieutordI 20021 : Passamonti et al. 2009al : Haskell et al. 20031 . The parameters for the two fluids are then 
given by En = XpS and Sp = e—Sn- For models A we must also specify the proton fraction and the sy mmetry energy term. Thes e 
are, respectively, set to Xp = 0.1 and a — 0. For a discussion of the effect of a on the spectrum, see lPassamonti et al.l (|2009al ). 
From the numerical simulations we determine the mode frequencies with an FFT of the time-evolved perturbation variables. 



In ord er to identify the differe n t mod es, we use also the eigenfunction extraction technique developed by IStergioulas et al, 
I 2004h and bimmelmeier et all \200di ). 

In Fig. Owe show, for the non-stratifled A models, some of the axisymmetric frequencies of the quasi-radial {I = 0) and 
quadrupole {I = 2) modes. In the left and right panels we show the "ordinary" and "superfluid" modes, respectively. These 
two mode families are decoupled in non-stratifled stars, and in fact the results in Fig. [5] do not hint at any interaction in 
the spectrum. However, within the sets of ordinary and superfluid modes avoiding crossings may appear. For instance, the 
ordinary quasi-radial mode H° and the ordinary pressure mode ^p§ seem to have an avoiding crossing when the star is rotating 
at 90% of the mass shedding limit. The effects of the chemical coupling on the spectrum is evident in Fig. [6] where we show 
some of the axisymmetric modes for the C models. In this case, the superfluid fundamental mode ^f ^ and the ordinary flrst 
pressure mode ^pi interact through an avoiding crossing near 90% of the Kepler limit. 

For a given I multipole, the non-axisymmetric modes of non-rotating stars have a degeneracy with respect to m. Rotation 
removes this degeneracy and splits each mode into 21 + 1 distinct branches. Besides the m = case considered above, 
we consider the |m| — 2 modes that have a pro- and retro-grade motion with respect to the star. In Fig. [71 we show the 
fr equencies of the acoustic m odes for models A (left panel) and C (right panel). The present results improve on the analysis 



of Passamonti et al 



(20093), that studied the dependence of the rotational splitting on the entrainment parameter within 
the Cowling approximation. The main improvement concerns the introduction of the gravitational potential perturbations. 
However, this does not alter the qualitative effects of rotation on the splitting. 
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Figure 8. This figure illustrates the efi'ect, in the weak drag regime, of the mutual friction on the stellar oscillations. For the model 
C2 with e = 0.5, the left panel displays the radial component of the variables 5Vnp (upper panel) and <5wnp (lower panel) for two long 
simulations with S = and B = 5 X 10~^, respectively. The horizontal axis shows the dimensionless evolution time. The lower- left panel 
shows that the counter-moving degrees of freedom are damped due to mutual friction dissipation. In the right panel, we show an FFT of 
the function (Si^^p for the K = 5 X 10~^ case. In order to study the mode amplitude variation with time, we have performed an FFT of 
the first part of the simulation, where C t(Gpo) ^ 500 (solid line), and subsequently of the second part, where 500 ^ t{Gpo) 1000 
(dashed line). On the horizontal axis is shown the dimensionless mode frequency a/ (Gpo)^^^, as measured in the rotating frame. As 
expected, the superfluid modes exhibit a faster damping than the ordinary modes. 



6.2 Mutual friction effects on the spectrum 



In order to study the effects of the mutual friction force on the oscillation spectrum, it is useful to write the momentum 
equation for the relative motion between protons and neutrons. We can combine the Euler-type equations (|32p and obtain 
the following expression in the rotating frame: 



(1 - e) dtSw'^" = -V (Sfip - Sp.n) - 2B'n x ^w^" + 2Z3 17 x n x Sw^"" 



where we have defined 



(78) 



(79) 



B' = 1 

Equation (|78|l makes the effects of the mutual friction parameters B and B' more evident. The term that includes B is 
dissipative and tends to damp the relative motion, and consequently mainly affects the superfluid modes. If the co- and 
counter-moving degrees of freedom are coupled, for instance due to the EoS, the mutual fri c tion dissipation affects also the 
ordinary mo des. Results to this effect, have been provided by (Lindblom fc Mendelll 1995 : Andersson et al. 20091 ) for the 
f-modes and jLindblom fc Mendelll [ioool : iLee fc Yoshidalbooj : iHaskell et al.ll2009l ) for the r-modes. The term proportional to 
B' modifies the Coriolis force, as one can see in equation (|78p . Its effects are not dissipative, but may change the frequencies 
of the superfluid modes. This is certainly expected in the case of the inertial modes as they are rotationally restored, but we 
will see that the non-axisymmetric fundamental modes can also be affected. 



20091 ) 



The magnitude of the mutual friction can be studied by introducing a dimensionless drag parameter TZ, deflned by ijHaskell et al 



B = 



7^ 



B' 



(80) 



1 + 7^2 ' 1 -f- 7^2 ■ 

Two extreme drag regimes can then be discerned for the mutual friction force. In the "weak" drag regime TZ <^ 1, whereas 
the "strong" drag regime corresponds to 7?. 3> 1. The most commonly considered cause of mutual friction is the scattering of 
the electrons off the magnetic field of the neutron vortices. This mechanism is firmly in the weak drag regime, where B <^ 1 
and B' B. In this case, we expect the mutual friction to act mainly on the mode damping. It should have negligible effects 
on the oscillation frequencies themselves. 



Recent discussions sugg est that the strong drag regime may lead to interesting, potentially important, results (jHaskell et al, 
20091 : I Andersson et al.ll2009l ). Since our level of theoretical understanding is not sufficient to rule out this case, we also consider 
the TZ ^ 1 regime. From equations (|80|) we see that in the strong drag regime S' ~ 1 and B <Si B' . The main effect should 
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Figure 9. This figure shows the effects of mutual friction, in the strong drag regime, on the rotational splitting of the superfluid 
I = m = 2 f-mode. The axis labels are shown in dimensionless units and the mode frequencies are determined in the rotating frame. For 
the sequence of models A with e = and Xp = 1/3, we study the B' = and B' = 1 cases. In the left panel, we show the effects of mutual 
friction on the pro-grade ^fp and retro-grade ^fj? modes, respectively, and the averaged frequencies (a) = (tTp -|- <7r) /2 of the two mode 
patterns. In the right panel, we show the deviation of the I = m = 2 mode defined by equation II83I I. The dashed lines are determined 
by using equation H84|l . These empirical relations agree very well with the values of the frequency deviation for the B' = I case, which 
are shown with filled circles. 



then be on the mode frequencies, while the dissipation can be considered negligible. In principle, we could explore also the 
intermediate regime, where TZ ~ 1 and both energy dissipation and frequency changes are important. However, this case is 
essentially a combination of the effects that we can study in the weak and strong regimes. Hence, we do not consider the 
intermediate regime in this work. 



6.2.1 Weak drag regime 

Let us first consider the weak drag regime by evolving in time the oscillations of model C2. In Fig. [S] we show the results from 
two long simulations where we have fixed TZ — and TZ = 5 x 10~^, respectively. In the left panel, we show the grid-averaged 
value of the velocities 5Vnp — 5vn + 5vp and 5wnp. In the upper- left panel, the two curves appear similar showing a weak 
damping that is mainly due to the numerical dissipation. In fact, the quantity iSVnp describes the evolution of the co-moving 
degree of freedom, which is weakly affected by the weak mutual friction. Looking more carefully at the results, we note that 
some damping is present in the TZ = 5 x 10"* case. This is due to the chemical coupling with the counter-moving degree of 
freedom, which is strongly damped. This is evident from the results in the lower-left panel of Fig. [H] which show that the 
amplitude of the relative velocity Swnp decreases during the evolution. 

These results suggest that, as expected, superfluid modes are damped faster than the ordinary modes. In order to study 
how the mode amplitude changes during the evolution, we divide the time-evolved data into two equal sets and perform an 
FFT for each part. Results for the variable (5wnp in the TZ — 5 x 10~* case are shown in the right panel of Fig. [8] We see that 
the superfluid fundamental and first pressure modes are dampe d faster then t heir o rdinary counterparts. 

The effect of the mutual friction has also been tested by ISiderv et al. l|2O10l) . by comparing the glitch spin-up time 



extracted by our numerical evolutions against an analytical formula derived within a body-averaged approximation. 

While our results demonstrate good progress, they are not quite satisfactory in one important respect. Ideally, one would 
like to be able to extract both oscillation frequency and damping time for the different modes seen in the evolution. However, 
so far we have not managed to extract the mutual friction damping rate of individual oscillation modes with the desired 
precision. This is basically because of the fact that the damping is very slow. It is also sensitive to the velocity lag between 
the two fluid components. At the present time it is not clear to us whether a time-evolution code provides a useful alternative 
to frequency-domain calculation for the damping-rate problem. 
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Figure 10. For model A2, we show (in the left panel) how the I = m = 2 superfluid f-mode depends on the proton fraction in the strong 
drag regime with B' = 1. The vertical axis displays the ratio of the quantity Act for the B' = 1 and B' = cases, see equation 1184 l l. In 
the right panel, we show the I = m = 2 ordinary ^r° and superfluid ■^r'^ modes for the sequence of models A with e = a = and TZ = 10^ 
{B' = 1). In the strong drag regime, the superfluid '^r'^ mode exhibits a strong dependence on the proton fraction Xp and thus on the 
parameter B' . 



6.2.2 Strong drag regime 

Next we explore the effects of the mutual friction in the strong drag regime, focussing on the I — m — 2 superfluid f- and 
r-modes. The parameter B' now dominates the mutual friction force aflFecting the Coriolis term in equation ((78}. 

The first aspect we want to understand is whether the rotational splitting of the superfluid f-mode is modified by the B' 
parameter. Based on our expectations, we assume that the frequency of the mode is described by the following relation 
up to order Q'^; 

a" ^cr'NR + ci{e,a,B',m)ft + 0{n^) , (81) 

where ci depends on the azimuthal index m and the stellar parameters e,a and B' . For B' = 0, we have already studie d 
the dependence of the ^f' mode on the entrainment parameter e and the symmetry energy term a ( Passamonti et al.ll2009al ). 



Therefore, we focus on the £ — a = case and vary the parameter B' . Using our previous results, we can re- write equation ifST)) 
as follows; 

cr' = a^fl + B'n + {if) , (82) 

We can then test this result against the numerical simulations. 

We first study a sequence of A models with a;p = 1/3, 7?. = 10'^, and B' — —2. According to equation (|82l) . we would 
expect the pro- and retro-grade mode-branches to be exchanged compared to the 7?. = case. In Fig. [9] we show the ^f° 
mode for stellar models A rotating up to the mass shedding limit with TZ — and TZ = 10^, respectively. The results show 
that (|82p describes the ^f' mode very well in the strong drag regime. We find that the scaling is quite accurate for stars up 
to il/Q,K ~ 0.9 (note that this analysis is not reported in Fig. [9]). However, the agreement is not so good when the mutual 
friction vanishes. It seems that for TZ = the effects of the centrifugal force becomes important for slower rotating models 
than in the TZ — 10'^ case. This behaviour is evident in FigO when we consider the averaged frequency between the m = 2 
pro- and retro-grade modes, i.e. (a) — ((jp -|- CTi) /2. 

However, we can take into account the effects of the centrifugal force on the average mode frequency and determine a 
connection between the superfluid f-mode frequencies in the strong and weak drag regimes. To this end, we define the mode 
deviation from its averaged value: 

Ao-= = 0-= - (a=) , (83) 
We then expect, from equation (I82p . that the following relation is valid: 

Aas'=i ~ B'Aas'^o ■ (84) 
In the left panel of Fig. |9l we show the quantity Aa for the ^f° mode in the strong drag regime and for vanishing mutual 
friction. The results for the TZ = 10^ case agree very well with the values obtained from equation (|84p . 

So far, we have studied a sequence of rotating stars with fixed proton fraction. Now, we test relation (|84|) by varying Xp 
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Table 3. This table provides the fitting coefficients 03,05 of equation II86I I. and their errors AC3 and AC5, for the superfluid ^r'' modes 
in the strong drag regime. The results correspond to the r-modes of four sequences of models A with TZ = 10^ and e = a = 0, and where 
the proton fraction takes the values shown in the first column. In the second column, we give the value of the parameter B' . For the first 
two models, we do not show the values of C5, as we fit the r-modc frequencies with C5 = 0. 





B' 


C3 


Ac3 


C5 

xl0~3 


Acs 
xlO-3 


2/5 


-1.5 


0.11338 


2.266 






1/3 


-2.0 


0.01121 


0.958 






1/4 


-3.0 


-0.05599 


0.948 


6.275 


0.233 


1/5 


-4.0 


-0.05982 


1.742 


4.013 


0.241 



and choosing the rotational model A2 with TZ = 
proton fraction is well described by equation (f84 



10' 



The results in Fig. \W[ show that the scaling of the ° mode with the 



However, when Xp = 0.1, there is a small difference between the numerical 

It is 



and analytical values. This effect might be due to the second order terms that we have neglected in the expansion (|82[1 
natural that these become important when the parameter B' is close to 10. 

Let us now study the behaviour of the superfluid I = m — 2 r-mode in the strong drag regime. Oscillations restored 
by t he Coriolis force generate t he class of inertial modes, which can be classified (by their parity) as axial-led or polar- 
led (|Lockitch &: Friedman|[l999h . The ordinary r-modes form a sub-set that is purely axial in the slow-rotation limit. The 
superfluid problem is somewhat different in that a purely axial superfluid r-mode exists only in non-stratified stars. When 
composition g radients are present , the superfluid r-mode acquires a polar component and assumes the nature of a general 



inertial mode I Haskell et al 



2009!). 



For a constan t density stellar mo del with B 
following relation ( Haskell et al. 20091 ) : 



■ 0, the frequency of the r'' mode in the rotating frame is described by the 



10/ o 



2m'yeB'9. 



l{l + l) 

where a° is the frequency of the ordinary r-mode. In the case of Z = m = 



(85) 

2 we have a° = 2Q/3. For compressible models, 
equation (|85|) approximately describes the frequency of the ^r^ mode only for slowly rotating stars. In fact, when a star 
rotates rapidly the effects of O (fi^) must be taken into account. In our time-evolutions, the rotational deformation of the 
star is completely described by the axisymmetric background. Meanwhile, in the slow-rotation approximation the equilibrium 
configuration remains spherical and the rotational effects on the spectrum are describe d by a perturbation expansion in Q. 
By using a slow-rotation approximation up to O (fi^) and the Cowling approximation, Haskell et al. ( 2009h determined the 
frequency of the '^r° and '^r^ modes in clo sed form. For the se quence of non-stratified A models with zero mutual friction, we 
have compared the r-mode frequencies of Haskell et al.l (I2OO9I ) with t he spectrum extracted by the time-evolutions and found 



by 1 

i). 



an agreement to better than 3% up to models with Q/Qk — 0.77 (jPassamonti et al.ll2009al ). For faster rotation, the slow- 
rotation approximation would require the calculation of terms of higher order than O which can be computationally 
prohibitive. In the strong drag regime, the effects of the higher order pertubative terms can become important even for 
relatively slowly rotating models, as a large value of B' increases the effective strength of the Coriolis force. 

We study the superfluid r-modes of the rotating A models, where we flx the values of the entrainment and symmetry energy 
to zero , £ = a — 0. The effects of these two parameters on the r-mode spectrum have already been studied by Passamonti et al.l 
1 2009al ). In this paper, we focus on the effects of the mutual friction parameter B' by choosing TZ = 10^ and consider four values 
of the proton fraction, namely Xp — 2/5, 1/3, 1/4, 1/5. The extraction of the r-mode frequencies from the time-evolutions 
requires longer simulations, as these modes are in the low-frequency regime. In order to save computational time, we adopt 
the Cowling approximation. In the right panel of Fig. [10] we show the ordinary ^r° and superfluid '^r'^ modes for different 
proton fractions. The ^1° mode has a retro-grade motion with respect to the star and is not affected by the parameter Xp. In 
contrast, the '^r'' mode depends strongly on Xp and has a pro-grade nature, as B' is negative. 

In order to understand the behaviour of the superfluid r-modes, we assume that for e = a = the frequency of a 
counter-moving I = m = 2 r-mode is described by the following relation; 



+ cJb' 



+ CS B' 



(86) 



where the frequencies and angular velocities are expressed in dimensionless units, while C3 and C5 are two fitting parameters. 
The numerical spectrum, shown in Fig. 1101 is well described by the first term of equation (|86p only for slowest rotating models. 
For stars with Xp > 1/4 the agreement is good up to — 0.4, while for a;p 1/4 the range reduces to Q/Qk ^ 0.2. In 

particular, we learn from the results in Fig. [lO] that the mode pattern changes concavity for increasing values of B' . This can 
be an effect of the O {^'^) and O (fi''') terms of equation (|86l) . Therefore, we fit our numerical data with equation (|86|) and 
determine the parameters C3 and C5. For models with > 1/4, a good fit can be determined by setting C5 = and calculating 
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Figure 11. This figure shows the waveform of the gravitational-wave signal for the C2 model with constant entrainment parameter 
e = 0.5. We show the evolution of the two independent initial conditions ID-N (dashed-line) and ID-P (solid-line). The axis labels are in 
dimensionless units. 



only the coefficient C3. For the ordinary ^r° mode we obtain C3 — 0.4852 ± 0.0126, while for the superfluid ^r'' modes the 
results are given in Table [S] When the proton fraction is smaller, i.e. Xp ^ 1/4, we must use the entire equation (|86p . The 
results of the corresponding fits are given in Table [S] In particular, when a;p ^ 1/4 we note a sign change in the parameter cs 
that represents the concavity variation of the mode pattern. 



6.3 Glitch gravitational signal 

We now turn to the gravitational-wave signal generated by initial axisymmetric configurations such that the protons and the 
neutrons rotate with a velocit y lag. As di scussed in Section [3. 2 1 these configurations can be determined with the perturbative 
approach developed bv lYoshida fc Eriguchi (2004). We have already considered this problem (Siderv et al. ,2010 ) in the context 
of pulsar glitch es. The following dis c ussion provides additional, more technical, details on these results. 



Within the lYoshida fc Eriguchil (|2004l ) approach all the initial axisymmetric, non-corotating configurations of a corotating 
background can be constructed as a linear combination of two independent classes of initial data, see Section [3.21 In the first 
class, only the neutrons move relative to the corotating background, i.e. (5f2n,5r2p) — (1,0), while in the second class only 
the proton velocity is different from the corotating background, {5fln, SQp) = (0, 1). We will refer to these two configurations 
as initial data N (ID-N) and P (ID-P), respectively. The initial mass density Sp^^, chemical potential Sjl^ and gravitational 
potential 5$ can be directly determined from equations (|25|l - (|26p for the two sets of initial data ID-N and ID-P. Meanwhile, 
for the velocity field perturbation we consider 

5vx = Sn^ {etc X r) . (87) 

These solutions can be rescaled to any required glitch size if we note that the crust spin-up can be associated with the proton 
velocity lag SQp. In fact, we expect an efficient coupling between the crust and the outer core protons due to the magnetic 
field, and we can then assume that the charged particles corotate. The rotational lag between superfluid neutrons and the 
protons can be estimated by considering angular momentum conservation: 

5J = InSQn + IpSnp + [Sin + Sin) fie = . (88) 

Here J is the total angular momentum, Jx is the moment of inertia of each fluid constituent, and its perturbation 5/x is 
defined by 



5L 



= / Sp^{r' sine fdr' . (89) 
Jo 



The initial relative velocity lag that describes a glitch is then given by 
^^^'^ - fie 



(90) 



Sfin = -^{ipSnp + sm,) , (91) 
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Figure 12. In this figure we illustrate the efl'ect of the Cowling approximation on the gravitational signal. The background star is the 
C2 model with e = 0.5, and the initial condition is that used for the glitch model, i.e. SQp = —10^'' and SQn = 7.74 X 10~*. The 
waveforms are shown in the left panel and the PSD in the right panel. The signals determined in Cowling approximation and with 
gravitational potential perturbation are shown with dashed and solid lines respectively. In the right panel we note the effect of the 
Cowling approximation on the acoustic modes. For the initial data considered in this work, the Cowling approximation model generates 
a larger gravitational- wave amplitude than in the case when the gravitational potential perturbation is accounted for. 



where SI = Sip + Sin- 

For each background model, we can evolve the two independent initial data sets ID-N and ID-P. If we consider a generic 
perturbation 5f for an arbitrary initial configuration, we can determine the evolution from the following linear combination: 

Sf = SfySnn + SfpQp , (92) 

where (5/n and Sfp are the perturbation variables related to ID-N and ID-P, respectively. 

We have evolved the ID-N and ID-P configurations for the C2 model with e — 0.5. In Fig. 1111 we show a part of the time 
evolution of the quantity A^" = Af-^ +Ap'^ determined from the stress- formula (|59[1 . Actually, we show the dimensionless quan- 
tity Al'^ = A'^^ / (^Gp^R^q) that is directly determined by the numerical code. For different values of the stellar parameters, the 
gravitational- wave amplitude can be calculated from equations (|71|) and (|73|) . The initial data ID-N generates a gravitational 
signal that is about an order of magnitude larger than the ID-P initial data. We have studied stellar models with different 
proton fraction and noticed that the amplitude difference between the ID-N and ID-P initial data scales with the proton 
fraction of the background model. This is expected as the dynamics of the mass constituents generates the gravitational-wave 
signal. 

For the glitch initial data, we study the effect of the Cowling approximation on the gravitational-wave signal. To do this, 
we consider two simulations for the same model C2 and entrainment parameter e = 0.5. The only difference is that, in one 
case we neglect the perturbation of the gravitational potential S^. In Fig. 1121 we show the time evolution of the quantity 
AI'^ and the related Power Spectrum Density (PSD), which is defined as PSD(A^'') — |A^"|. In the Cowling approximation, 
we extract the gravitational signal with the momentum formula (|58p . as the stress formula (|59p is not well defined when 
(5$ = 0. From the results in the left panel of Fig. 1121 we note that the Cowling approximation generates a signal that is about 
five times larger than the result when 5$ is included. Furthermore, as expected, the Cowling approximation introduces a 
deviation in the mode frequencies. This difference is evident in the right panel of Fig. 1121 where the error is about 38% for 
the fundamental quasi-radial mode (F°), 32% for the axisymmetric / — 2 f-mode (^f°), and 10% for the first pressure mod e 
(^Pi). These results agree well with the results of similar comparisons ( Yoshida fc Koiima 1997 : Yoshida fc Eriguchil 2001 ). 
Regarding the amplitude of the gravitational-wave signal, we find that the relative oscillation amplitude between the Cowling 
approximation and the full problem depends on the initial data. Hence, this result is not generic. 



7 CONCLUSIONS AND DISCUSSION 

We have studied the dynamics of superfluid rotating neutron stars, focussing on the nature of the oscillation spectrum, the 
effects of the mutual friction force on the oscillations and the hydrodynamic spin-up phase of pulsar "glitches" . Adopting the 
Newtonian two-fluid model, we evolved in time the perturbed dynamical equations on axisymmetric equilibrium configurations. 
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This approach allows us to derive the spectrum of axisymmetric and non-axisymmetric oscillation modes of stellar models 
that rotate up to the mass shedding limit. In this work, we have improved on previous studies by including the gravitational 
perturbation and the mutu al friction force. The spe ctrum is then determined with a better accuracy, as we no longer use 
the Cowling approximation ( Passamonti et al. 2009ah . From the computational point of view, we have to solve the perturbed 
Poisson equation together with the linearised momentum and mass conservation equations. We have numerically evolved the 
hyperbolic equations with a Mac-Cormack algorithm, while the elliptic equation for the gravitational potential is solved at 
each time step with a spectral method. 

In our current model the rotating background models are pure fluid, i.e. without an elastic crust region, neutrons and 
protons corotate and are in /3-equilibrium. In superfluid stars, the co- and counter-phase motion of the two fluid constituents 
can be coupled by composition gradients and this influences the dynamics. In order to consider this effect we have studied 
two simple polytropic equations of state that generate distinct sequences of stratified and non-stratified rotating stars. These 
background models are simplistic, and we must improve on this aspect if we want to decode the complexity of astrophysical 
observations. Certainly, we must add an elastic crust to the model and relax the co-rotation assumption between the two 
fluids. If we want to use more realistic equations of state we also need to translate the model to General Relativity. We are 
currently working on all these issues. 

In neutron stars, the mutual friction force may have both dissipative and non-dissipative effects. The dissipative part of 
the force, which is dominant in the weak drag regime, mainly damps an oscillation mode. Meanwhile, the non-dissipative term 
dominates in the strong drag regime, essentially modifying the oscillation spectrum. We have studied the two drag regimes 
and showed that our numerical code effectively reproduces the mutual friction damping of the two-fluid relative motion. For 
non-stratified stars, the co- and counter-moving degrees of freedom are uncoupled and only the superfluid modes are damped. 
When the stellar m odel is stratif i ed, th e damping affects also the ordinary modes. The accuracy of our numerical code has 
also been tested in ISiderv et al. where determined the glitch sp n-up time and compared it to a simple analytic 

formula. However, we are not yet able to extract (with useful precision) the mutual friction damping time of individual 
oscillation modes from our numerical evolutions. More work is needed to establish to what extent one should expect to do 
this within our computational framework. For the strong drag regime, we have studied the effect of the mutual friction and 
composition variation on the rotational splitting of the superfluid I — m = 2 f-mode and on the frequencies of the I = m — 2 
superfluid r-mode. The main effect is a change of propagation direction of the modes with respect to the background rotation. 
A mode that is pro-grade (retro-grade) in the weak drag regime may become retro-grade (pro-grade) in the strong drag 
regime. We have determined the numerical frequencies of the f- and r-modes for the rotating sequence of non-stratified stellar 
models and provided simple empirical expressions based on the numerical data. For constant mutual friction parameters, the 
non-axisymmetric splitting of the superfluid f-mode and the r-mode frequencies depends on the inverse of the proton fraction. 



Finally, we provided relevant technical details for the hydrodynamical models for pulsar glitches discussed bv lSiderv et al 



The initial conditions for the glitch evolutions describe two fluids that rota te with a small veloc i ty lag . These config- 
urations were been determined using a perturbative approach first introduced by Yoshida fc Eriguchil ( 2004h . We extended 
this method to implement different EoS and consider non-corotating initial configurations that conserve the mass of each 
fiuid constituent. Moreover, we derived the detailed quadrupole gravitational extraction formulae for / = 2 oscillation modes 
of a superfiuid star. We determined the perturbative expressions for the momentum and stress formulae that can be used to 
improve the numerical extraction of the gravitational-wave signal (reducing the order of the time derivative of the standard 
quadrupole formula). Wc determined the gr avitat ional- wave strain for the two independent initial glitch configurations that 
are obtained with the lYoshida fc Eriguchil (j2004h approach. For a given background rotation, these results can be used to 
estimate the gravitational signal for any glitch size. Furthermore, we have showed the effect of the Cowling approximation on 
the glitch gravitational-wave strain and the oscillation spectrum. 

With the progress described in this paper, our programme of studying superfluid neutron star dynamics by time-evolutions 
of the linearised equations has reached the point where we need to add key physics to the model. The natural step would be 
to account for the elastic neutron star crust with the expected interpenetrating neutron superfiuid. This requires us to change 
the computational framework somewhat, as it is natural to discuss the elasticity in term of Lagrangian perturbation theory. 
Moreover, we need to address various issues associated with vortex pinning by the crust nuclei. This problem requires additional 
force contributions at the level of individual vortices, and we need to develop a suitable smooth-averaged hydrodynamics 
description if we want to make progress. We are currently working on both these issues. It would also be relevant to extend 
our models to general relativity. This is essential if we want to be able to use realistic supranuclear equations of state. As 
long as we make use of the relativistic analogue of the Cowling approximation this generalisation should be straightforward, 
but if we want to account for the dynamics of spacetime the problem becomes much more involved. If we want to consider 
realistically "layered" neutron stars we also need to improve our understanding of the different phase-transitions, e.g. in the 
vicinity of the critical density/temperature for the onset of superfluidity, and how these regions affect the large scale dynamics. 
We face a number of challenging questions, but there is no reason why we should not be able to resolve the relevant issues 
and progress towards the construction of realistic dynamical neutron star models. 
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APPENDIX A: GW EXTRACTION 

In this Appendix we determine the momentum formula (|65p and the stress formula (|66|l for the {I, m) = (2, 2) gravitational 
signal. For the axisymmetric {I, m) — (2, 0) case, we have used the perturbative version of the momentum and stress formulae 
used bv lFinn fc Evans! !\im(ji ). 

The aim is to reduce the order of the time derivatives in the quadrupole gravitational-wave formula. To do this we consider 
the quantity (|63[l : 

A^^ = -^ I drSpr^Y;2. (Al) 
With the use of the mass conservation equations of each fluid component: 

ftpx + V, (pxt;;) = , (A2) 



we can determine the momentum-formula where only a first order time derivative appears. When we perturb equation (|A2 
and introduce it in (IA1|) we obtain: 



(A3) 



where in the last step we have used the Gauss Theorem. After some calculation, equation (IA3|I leads to the expression (|65|l . 

With a similar method, we can determine the stress-formula and eliminate the time derivatives from the quadrupole 
formula. In this case, we must use the momentum conservation equation that for a superfluid component is given by 

dt {p^pD + Vfe (^pWy,Pi^ + /5xVi/ix + PxV,* -I- pxExiffe^ViUx = , (A4) 

where the momentum of the fluid component is deflned as follows: 

P? = Vi + Exwf . (A5) 

For a two-fluid model with neutron and proton as components, the total momentum equation is then given by the following 
expression: 

where we have used the definition of the the generalized pressure l|Prixll2004) : 



(A6) 



V* = PnV/in + PpV/ip - -pxExV (lUpn) , (A7) 



and we have re-written the gravitational potential term by using the Poisson equation: 



pV,$ = ^V'' (^Vfc$V.$ - ^V,<E>V^<E>j . (A8) 
Perturbing equation (|A6|) and considering a corotating equilibrium configuration, i.e. lUnp = 0, we obtain: 

^ 5 (pnv; + ppvf) = -VkS (pn^Sp? + Pp«Jpr) - - -^^"S (^Vfe-I>V,$ - ^ v,$V^$^ , (A9) 
where now for corotating background the pressure perturbation is given by 

VSP = S (pnVpn + PpVpp) . (AlO) 

We can now introduce equation (|A9|) into equation (|A3I) and use the Gauss theorem. We obtain: 

^ fdrSp = fdrS fpnvtpn + pAPp + i^'*^' ^2*2) (All) 



where both the pressure and the last term of equation (|A9|) vanish, as V (r Y22) = 0. After some further calculation, we can 
derive equation (|66p from HA11[I . 



24 A. Passamonti & N. Andersson 
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